All the MED nodes
physeq <- ps.filter
colnames(tax_table(physeq))[7] <- "Strain"
levels(physeq@sam_data$Species) <- c("AA", "CP", "CQ")
# levels(physeq@sam_data$Species)
levels(physeq@sam_data$Strain) <- c("B", "CE", "G", "L", "W-")
# levels(physeq@sam_data$Strain)
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Genus",
tax_add = NULL,
ntax = 100) -> p
## Loading required package: ampvis2
## Warning: There are only 13 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free")+
scale_fill_viridis_c(breaks = c(0, 0.01, 0.1, 1, 10, 100),
labels = c(0, 0.01, 0.1, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma=0.001), # Scaling factor for the linear part of pseudo-log transformation
na.value = 'transparent') -> p1
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p1[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p1
Wolbachia nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Wolbachia") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 20 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p2
Asaia nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Asaia") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 24 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p3
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p3[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p3
Elizabethkingia nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Elizabethkingia") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 5 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p4
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p4[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p4
Legionella nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Legionella") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p5
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p5[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p5
Chryseobacterium nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Chryseobacterium") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 4 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p6
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p6[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p6
Erwinia nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Erwinia") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p7
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p7[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p7
Serratia nodes
physeq %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus == "Serratia") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Strain", "Organ", "Field", "Species" ),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
## Warning: There are only 2 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p8
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p8[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p8
Check Blanks
# import phyloseq object before decontam
setwd(path_rdata)
ps <- readRDS("MED_phyloseq.rds")
colnames(tax_table(ps))[7] <- "Strain"
# 70 nodes (all)
ps %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Control"),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p
p + facet_grid( ~ Control, scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p9
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p9[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p9
# Pseudomonas nodes (3)
ps %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus=="Pseudomonas") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Control"),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p10
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p10[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p10
# Enhydrobacter nodes (3)
ps %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus=="Enhydrobacter") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Control"),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p11
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p11[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p11
# Rahnella1 nodes (3)
ps %>%
transform_sample_counts(function(x) x/sum(x) *100) %>%
subset_taxa(Genus=="Rahnella1") %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Control"),
tax_aggregate = "Species",
tax_add = NULL,
ntax = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
theme(axis.text.y = element_text(face="italic", angle = 0, size = 12))+
scale_fill_viridis_c(breaks = c(0, 0.01, 1, 10, 100),
labels = c(0, 0.01, 1, 10, 100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p12
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p12[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p12
Save heatmaps
setwd(path_plot)
tiff("1Ed_MED_heatmap.tiff", units="in", width=28, height=16, res=300)
p1
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_wolbachia.tiff", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_asaia.tiff", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_elizabethkingia.tiff", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_legionella.tiff", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_chryseobacterium.tiff", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_erwinia.tiff", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_serratia.tiff", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_blanks_all.tiff", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_blanks_pseudomonas.tiff", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_blanks_enhydrobacter.tiff", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen
## 2
tiff("1Ed_MED_heatmap_blanks_rahnella1.tiff", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap.png", units="in", width=22, height=10, res=300)
p1
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_wolbachia.png", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_asaia.png", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_elizabethkingia.png", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_legionella.png", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_chryseobacterium.png", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_erwinia.png", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_serratia.png", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_blanks_all.png", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_blanks_pseudomonas.png", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_blanks_enhydrobacter.png", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen
## 2
png("1Ed_MED_heatmap_blanks_rahnella1.png", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen
## 2